############################################
# Replication Script for Home and Away: Explaining the Paradoxical Political Attitudes of Indian Americans
# Authors: Sumitra Badrinathan with Annabel Richter and Josiah Gottfried
# Date: February 13, 2026
############################################

##### SETUP #####

# setwd() # set your working directory

library(pacman)

# loading in libraries
p_load(scales, ggplot2, ggpubr, readxl, grid, gridExtra, dplyr, ggthemes, srvyr, survey, 
       tidyverse, stargazer, pacman, install = TRUE)

# loading survey data
iaas <- read_csv("JHOP0001_OUTPUT_order.csv")
iaas_wtd <- as_survey_design(iaas, weights = weight)

# load states data
fec_states <- read_xlsx("FEC_2020.xlsx")
fec_states <- fec_states%>% 
  mutate(party = case_when(`Biden_D` > `Trump_R` ~ "D", 
                           `Trump_R` > `Biden_D` ~ "R"), 
         lean = case_when(party == "R" ~ "Red state", 
                          TRUE ~ "Blue state"))

##### CREATE VARIABLES #####

# creating variables for use later 
iaas_wtd <- iaas_wtd %>%
  mutate(
    # party variables
    usparty_ = factor(case_when(pid3 == "Other" ~ NA_character_,
                                TRUE ~ pid3), 
                      levels = c("Democrat", "Independent", "Republican", "Not sure")),
    inparty_ = factor(case_when(B5_Q3 == "Bharatiya Janata Party (BJP)" ~ "BJP",
                                B5_Q3 == "Don't know" ~ "Don't know", 
                                TRUE ~ "Not BJP"),
                      levels = c("BJP", "Not BJP", "Don't know")),
    modi_ = B5_Q4,
    modidem_ = factor(case_when(modi_ == "Approve" & usparty_ == "Democrat" ~ "Modi Democrat",
                                modi_ == "Approve" ~ "Other Modi",
                                usparty_ == "Democrat" ~ "Other Democrat",
                                TRUE ~ "Neither"),
                      levels = c("Modi Democrat", "Other Modi", "Other Democrat", "Neither")),
    modi_democrat = case_when(modi_ == "Approve" & usparty_ == "Democrat" ~ 1, 
                              TRUE ~ 0),
    
    # demographic variables
    religion_ = factor(case_when(religpew == "Hindu" ~ "Hindu",
                                 religpew == "Muslim" ~ "Muslim",
                                 TRUE ~ "Other")),
    whitenatthreat_ = factor(case_when(str_detect(B4US_Q23, " agree") ~ "Agree",
                                       TRUE ~ "Disagree"), 
                             levels = c("Disagree", "Agree")),
    hindunatthreat_ = factor(case_when(str_detect(B5_Q16, " agree") ~ "Agree", 
                                       TRUE ~ "Disagree"), 
                             levels = c("Disagree", "Agree")),
    education_ = factor(case_when(!str_detect(educ, "4-year|Post-grad") ~ "<Bachelor's",
                                  TRUE ~ educ)),
    immigrant_ = factor(case_when(B1_Q2 == "Yes" ~ "Born in USA",
                                  TRUE ~ "Immigrant")),
    cohort_ = factor(case_when(as.integer(B1_Q5) >= 2014 ~ "2014 or later",
                               as.integer(B1_Q5) < 2004 ~ "Before 2014",
                               TRUE ~ NA_character_)),
    state_lean = fec_states$lean[match(inputstate, fec_states$State)],
    state_ = factor(state_lean, levels = c("Red state", "Blue state")),
    income_ = factor(case_when(str_detect(faminc_new, "[0-9]{3},[0-9]{3} ") ~ ">$100,000",
                               faminc_new == "Prefer not to say" ~ NA_character_,
                               TRUE ~ "<$100,000")),
    gender_ = factor(case_when(gender == "Male" ~ "Male",
                               gender == "Female" ~ "Female",
                               TRUE ~ NA_character_)),
    age_ = factor(case_when(2020 - birthyr < 30 ~ "<30",
                            2020 - birthyr < 50 ~ "30-50",
                            TRUE ~ ">50"),
                  levels = c("<30", "30-50", ">50")),
    engagement_ = rowSums(select(.[["variables"]], B4Filler_Q4_1:B4Filler_Q4_4) %>%
                            mutate(across(everything(), ~if_else(.x == "selected", 1, 0)))),
    engagement_ := factor(case_when(engagement_ == 0 ~ "Low engagement",
                                    engagement_ == 1 ~ "Medium engagement",
                                    TRUE ~ "High engagement"),
                          levels = c("Low engagement", "Medium engagement", "High engagement")),
    identity_ = factor(case_when(str_detect(B4Filler_Q8, "unimportant") ~ "Unimportant",
                                 TRUE ~ B4Filler_Q8),
                       levels = c("Unimportant", "Somewhat important", "Very important")))

# dependent variable (difference context) 
iaas_wtd <- iaas_wtd %>%
  # B4US_Q18: cracking down on government critics (US), Q19: stricter immigration (US), Q21: Muslim travel ban (US), Q22: police force during peaceful protests (US)
  # B5_Q11: cracking down on government critics (IN), Q12: pro-NRC to deport immigrants (IN), Q13: pro-CAA to prevent Muslim immigrants (IN), Q15: police force during peaceful protests (IN)
  mutate(across(c(B4US_Q18:B4US_Q19, B4US_Q21:B4US_Q22, B5_Q11:B5_Q13, B5_Q15), list(code = ~case_when(.x == "Strongly support" ~ 1,
                                                                                                       .x == "Somewhat support" ~ 2,
                                                                                                       .x == "Somewhat oppose" ~ 3,
                                                                                                       .x == "Strongly oppose" ~ 4))),
         # B4US_Q20: affirmative action (US), Q14: affirmative action (IN) 
         across(c(B4US_Q20, B5_Q14), list(code = ~case_when(.x == "Strongly oppose" ~ 1,
                                                            .x == "Somewhat oppose" ~ 2,
                                                            .x == "Somewhat support" ~ 3,
                                                            .x == "Strongly support" ~ 4))),
         across(c(B5_Q11_code:B5_Q15_code, B5_Q14_code), ~-.x)) %>%
  mutate(diff = rowSums(select(.[["variables"]], c(B4US_Q19_code, B4US_Q21_code, B5_Q12_code, B5_Q13_code))),
         diff2 = rowSums(select(.[["variables"]], ends_with("_code"))))

fct_levels <- c("Treating members of all religious groups equally",
                "More permissive policies toward undocumented immigrants",
                "Preventing police use of force against peaceful protesters",
                "Protection of media from government censorship",
                "Affirmative action in university admissions")

##### FIGURE BY FIGURE RESULTS #####

##### figure 2 #####

  fig2_df <- iaas_wtd %>%
  mutate(modi_ = str_replace(modi_, "([Aa]pprove$)", "\\1\nof Modi") %>%
           as_factor) %>%
  group_by(interact(usparty_, modi_)) %>%
  summarize(prop = survey_prop())

h_total <- fig2_df %>%
  group_by(usparty_) %>%
  summarize(prop = sum(prop)) %>%
  mutate(modi_ = factor("Total", levels = c(levels(fig2_df$modi_), "Total")))

v_total <- fig2_df %>%
  group_by(modi_) %>%
  summarize(prop = sum(prop)) %>%
  mutate(usparty_ = factor("Total", levels = c(levels(fig2_df$usparty_), "Total")))

ggplot(data = fig2_df, aes(usparty_, modi_, label = percent(prop, accuracy = 1))) +
  geom_tile(aes(fill = prop), color = "white") +
  geom_label(aes(label = percent(prop, accuracy = 1)), size = 5.5) +
  geom_point(data = h_total, aes(color = prop), size = 30, shape = 19) +
  geom_label(data = h_total, aes(label = percent(prop, accuracy = 1)), size = 5.5) +
  geom_point(data = v_total, aes(color = prop), size = 30, shape = 19) +
  geom_label(data = v_total, aes(label = percent(prop, accuracy = 1)), size = 5.5) +
  scale_x_discrete(position = "top", limits = levels(v_total$usparty_)) +
  scale_y_discrete(limits = rev(levels(h_total$modi_))) +
  scale_fill_distiller(palette = "Blues", direction = 1) +
  scale_color_distiller(palette = "Blues", direction = 1) +
  labs(caption = "<1% of respondents identify with another U.S. party") +
  theme_void() %+replace%
  theme(title = element_text(size = 16, hjust = 0),
        plot.caption = element_text(size = 12, hjust = 1, vjust = 1),
        axis.text = element_text(size = 14),
        legend.position = "none")


##### figure 3 #####
noother <- subset(iaas_wtd, iaas_wtd$variables$religion_ == "Hindu" | iaas_wtd$variables$religion_ == "Muslim")

noother$variables %>%
  ggplot(aes(B4US_Q24N_2, B5_Q18N_4, color = religion_)) +
  geom_point(alpha = 0.5, position = "jitter", size = 1.8) +
  scale_color_manual(values = c("Hindu" = "#F8766D", "Muslim" = "forestgreen")) + 
  xlab("Feelings toward Democratic party") +
  ylab("Feelings toward Modi") +
  theme_minimal() %+replace%
  theme(plot.title = element_text(size = 16, hjust = 0),
        axis.title = element_text(size = 14),
        axis.text = element_text(size = 12),
        legend.title = element_blank(),
        legend.text = element_text(size = 12),
        panel.grid = element_blank())



##### table 2 #####
india_rbd <- iaas_wtd %>%
  select(religpew, B5_Q11:B5_Q15) %>%
  mutate(across(-c(religpew, B5_Q14), ~ case_when(str_detect(., "oppose") ~ 1,
                                                  TRUE ~ 0)),
         B5_Q14 := case_when(str_detect(B5_Q14, "support") ~ 1,
                             TRUE ~ 0)) %>%
  group_by(religpew) %>%
  summarize(across(everything(), survey_mean, vartype = "ci")) %>%
  filter(religpew %in% c("Hindu", "Muslim")) %>%
  pivot_longer(cols = -religpew, names_pattern = "(.{6})(.*)", names_to = c("variable", "limit"), values_to = "value") %>%
  mutate(se := case_when(limit == "" ~ "prop",
                         limit == "_low" ~ "lower",
                         TRUE ~ "upper"),
         variable := factor(case_when(variable == "B5_Q13" ~ "Treating members of all religious groups equally",
                                      variable == "B5_Q12" ~ "More permissive policies toward undocumented immigrants",
                                      variable == "B5_Q15" ~ "Preventing police use of force against peaceful protesters",
                                      variable == "B5_Q11" ~ "Protection of media from government censorship",
                                      variable == "B5_Q14" ~ "Affirmative action in university admissions"),
                            levels = fct_levels)) %>%
  pivot_wider(id_cols = c(religpew, variable), names_from = se, values_from = value) %>%
  arrange(variable)

full_join(filter(india_rbd, religpew == "Hindu"), 
          filter(india_rbd, religpew == "Muslim"), 
          by = "variable", suffix = c("_hindu", "_muslim")) %>%
  mutate(difference = prop_hindu - prop_muslim,
         significant = case_when(lower_hindu > upper_muslim | upper_hindu < lower_muslim ~ "Yes",
                                 TRUE ~ "No")) %>%
  select(variable, significant, prop_hindu, prop_muslim, difference)

##### table 3 #####
usa_rbd <- iaas_wtd %>%
  select(religpew, B4US_Q18:B4US_Q22) %>%
  mutate(across(-c(religpew, B4US_Q20), ~ case_when(str_detect(., "oppose") ~ 1,
                                                    TRUE ~ 0)),
         B4US_Q20 := case_when(str_detect(B4US_Q20, "support") ~ 1,
                               TRUE ~ 0)) %>%
  group_by(religpew) %>%
  summarize(across(everything(), survey_mean, vartype = "ci")) %>%
  filter(religpew %in% c("Hindu", "Muslim")) %>%
  pivot_longer(cols = -religpew, names_pattern = "(.{8})(.*)", names_to = c("variable", "limit"), values_to = "value") %>%
  mutate(se := case_when(limit == "" ~ "prop",
                         limit == "_low" ~ "lower",
                         TRUE ~ "upper"),
         variable := factor(case_when(variable == "B4US_Q21" ~ "Treating members of all religious groups equally",
                                      variable == "B4US_Q19" ~ "More permissive policies toward undocumented immigrants",
                                      variable == "B4US_Q22" ~ "Preventing police use of force against peaceful protesters",
                                      variable == "B4US_Q18" ~ "Protection of media from government censorship",
                                      variable == "B4US_Q20" ~ "Affirmative action in university admissions"),
                            levels = fct_levels)) %>%
  pivot_wider(id_cols = c(religpew, variable), names_from = se, values_from = value) %>%
  arrange(variable)

full_join(filter(usa_rbd, religpew == "Hindu"), 
          filter(usa_rbd, religpew == "Muslim"), 
          by = "variable", suffix = c("_hindu", "_muslim")) %>%
  mutate(difference = prop_hindu - prop_muslim,
         significant = case_when(lower_hindu > upper_muslim | upper_hindu < lower_muslim ~ "Yes",
                                 TRUE ~ "No")) %>%
select(variable, significant, prop_hindu, prop_muslim, difference)


##### figure 4 #####
iaas_wtd$variables %>%
  filter(religpew %in% c("Hindu", "Muslim")) %>%
  ggplot(aes(x = diff2, y = after_stat(density), fill = religpew)) +
  geom_histogram(binwidth = 1, alpha = 0.4, position = "identity") +
  ggtitle("Distribution of difference scores for Hindus and Muslims", subtitle = "Difference score calculated from all questions") +
  scale_x_continuous(name = "Difference (higher number signifies more liberal in US context)") +
  scale_y_continuous(labels = percent) +
  theme_minimal() %+replace%
  theme(plot.title = element_text(size = 16, hjust = 0),
        plot.subtitle = element_text(size = 12, hjust = 0),
        plot.title.position = "plot",
        axis.title.x = element_text(size = 14),
        axis.title.y = element_blank(),
        axis.text = element_text(size = 12),
        legend.title = element_blank(),
        legend.text = element_text(size = 12))

##### table 4 #####
# redefining IVs
iaas_wtd <- iaas_wtd %>%
  mutate(
    # demographic variables
    religion_ = factor(case_when(religpew == "Hindu" ~ "Hindu",
                                 religpew == "Muslim" ~ "Muslim",
                                 TRUE ~ "Other"), 
                       levels = c("Muslim", "Hindu", "Other")),
    education_ = factor(case_when(!str_detect(educ, "4-year|Post-grad") ~ "<Bachelor's",
                                  TRUE ~ educ)),
    immigrant_ = factor(case_when(B1_Q2 == "Yes" ~ "Born in USA",
                                  TRUE ~ "Immigrant")),
    state_lean = fec_states$lean[match(inputstate, fec_states$State)],
    state_ = factor(state_lean, levels = c("Red state", "Blue state")),
    income_ = factor(case_when(str_detect(faminc_new, "[0-9]{3},[0-9]{3} ") ~ ">$100,000",
                               faminc_new == "Prefer not to say" ~ NA_character_,
                               TRUE ~ "<$100,000")),
    gender_ = factor(case_when(gender == "Male" ~ "Male",
                               TRUE ~ "Not Male"), 
                     levels = c("Not Male", "Male")),
    age_ = factor(case_when(2020 - birthyr < 30 ~ "<30",
                            2020 - birthyr < 50 ~ "30-50",
                            TRUE ~ ">50"),
                  levels = c("<30", "30-50", ">50")), 
    biden2020_ = factor(case_when(B2_Q4 == "Joe Biden" ~ "Voted for Biden", 
                                  TRUE ~ "Did not vote for Biden"), 
                        levels = c("Did not vote for Biden", "Voted for Biden")))

# creating DV 
iaas_wtd <- iaas_wtd %>%
  mutate(across(c(B4US_Q18:B4US_Q19, B4US_Q21:B4US_Q22, B5_Q11:B5_Q13, B5_Q15), list(code = ~case_when(.x == "Strongly support" ~ 1, 
                                                                                                       .x == "Somewhat support" ~ 2,
                                                                                                       .x == "Somewhat oppose" ~ 3,
                                                                                                       .x == "Strongly oppose" ~ 4))),
         across(c(B4US_Q20, B5_Q14), list(code = ~case_when(.x == "Strongly oppose" ~ 1,
                                                            .x == "Somewhat oppose" ~ 2,
                                                            .x == "Somewhat support" ~ 3,
                                                            .x == "Strongly support" ~ 4))),
         across(c(B5_Q11_code:B5_Q15_code, B5_Q14_code), ~-.x)) %>%
  mutate(diff = rowSums(select(.[["variables"]], c(B4US_Q19_code, B4US_Q21_code, B5_Q12_code, B5_Q13_code))),
         diff2 = rowSums(select(.[["variables"]], ends_with("_code"))))

redef <- svyglm(diff2 ~ religion_ + education_ + immigrant_ + state_ + income_ + gender_ + age_ + biden2020_, 
                design=iaas_wtd)

stargazer(redef, type = "text", 
          covariate.labels = c("Religion: Hindu", "Religion: Other", "Education: Undergrad", "Education: Post-grad", "Immigrant", "Blue state", "Income over 100k", "Male", "Age 30 to 50", 
                               "Age over 50", "Biden voter (2020)"), 
          dep.var.labels = "Context Difference", star.cutoffs = c(0.05, 0.01, 0.001))


##### figure 5 #####
nationalism <- iaas_wtd %>%
  select(religpew, B4US_Q23, B5_Q16) %>%
  mutate(across(B4US_Q23:B5_Q16, ~case_when(str_detect(.x, " agree") ~ 1,
                                            TRUE ~ 0))) %>% 
  group_by(religpew) %>%
  summarize(`White Supremacy-` = survey_mean(B4US_Q23, vartype = "ci"), `Hindu Majoritarianism-` = survey_mean(B5_Q16, vartype = "ci")) %>%
  filter(religpew %in% c("Hindu", "Muslim")) %>%
  pivot_longer(cols = -religpew, names_pattern = "(.*)-(.*)", names_to = c("variable", "limit"), values_to = "value") %>%
  mutate(limit := case_when(limit == "" ~ "prop",
                            limit == "_low" ~ "lower",
                            TRUE ~ "upper")) %>%
  pivot_wider(id_cols = c(religpew, variable), names_from = limit, values_from = value)

# plot
nationalism %>%
  ggplot(aes(x = religpew, y = prop, label = paste0(percent(prop, accuracy = 1)), color = religpew)) + 
  geom_point(size = 4) +
  geom_errorbar(aes(ymin = lower, ymax = upper), linewidth = 0.75, width = 0) +
  geom_text(nudge_x = 0.35, size = 5, show.legend = F) +
  facet_wrap(~variable, ncol = 2) +
  ggtitle("Percent who agree that each form of nationalism is a threat to minorities\n") +
  scale_y_continuous(labels = percent) +
  scale_color_manual(values = c("#F8766D", "#619CFF")) +
  theme_bw() %+replace%
  theme(plot.title = element_text(size = 14, hjust = 0),
        legend.title = element_blank(),
        legend.text = element_text(size = 12),
        strip.text = element_text(size = 12, margin = margin(8), vjust = 2),
        axis.title = element_blank(),
        axis.text.x = element_blank(),
        axis.text.y = element_text(size = 12),
        axis.ticks.x = element_blank(),
        panel.grid.minor = element_blank())


##### figure 6 #####

fct_levels <- c("Treating members of all religious groups equally",
                "More permissive policies toward undocumented immigrants",
                "Preventing police use of force against peaceful protesters",
                "Protection of media from government censorship",
                "Affirmative action in university admissions")

generic_rbd <- iaas_wtd %>%
  select(religpew, B3_Q1:B3_Q5) %>%
  mutate(across(B3_Q1:B3_Q4, ~ case_when(str_detect(., " important|support") ~ 1,
                                         TRUE ~ 0)),
         B3_Q5 := case_when(str_detect(B3_Q5, "oppose") ~ 1,
                            TRUE ~ 0)) %>%
  group_by(religpew) %>%
  summarize(across(everything(), survey_mean, vartype = "ci")) %>%
  filter(religpew %in% c("Hindu", "Muslim")) %>%
  pivot_longer(cols = -religpew, names_pattern = "(.{5})(.*)", names_to = c("variable", "limit"), values_to = "value") %>%
  mutate(se := case_when(limit == "" ~ "prop",
                         limit == "_low" ~ "lower",
                         TRUE ~ "upper"),
         variable := factor(case_when(variable == "B3_Q4" ~ "Treating members of all religious groups equally",
                                      variable == "B3_Q2" ~ "More permissive policies toward undocumented immigrants",
                                      variable == "B3_Q5" ~ "Preventing police use of force against peaceful protesters",
                                      variable == "B3_Q1" ~ "Protection of media from government censorship",
                                      variable == "B3_Q3" ~ "Affirmative action in university admissions"),
                            levels = fct_levels)) %>%
  pivot_wider(id_cols = c(religpew, variable), names_from = se, values_from = value) %>%
  arrange(variable)

full_join(filter(generic_rbd, religpew == "Hindu"), 
          filter(generic_rbd, religpew == "Muslim"), 
          by = "variable", suffix = c("_hindu", "_muslim")) %>%
  mutate(difference = prop_hindu - prop_muslim,
         significant = case_when(lower_hindu > upper_muslim | upper_hindu < lower_muslim ~ "Yes",
                                 TRUE ~ "No")) %>%
  select(variable, significant, prop_hindu, prop_muslim, difference)


all_rbd <- bind_rows(
  generic_rbd %>% mutate(context = "generic"),
  usa_rbd %>% mutate(context = "usa"),
  india_rbd %>% mutate(context = "india")
)

all_rbd %>%
  # modify variables so that labels look good on the chart
  group_by(variable, context) %>%
  mutate(greater = (prop == max(prop)),
         context = factor(case_when(context == "generic" ~ "Generic",
                                    context == "usa" ~ "USA",
                                    context == "india" ~ "India"),
                          levels = c("Generic", "USA", "India")),
         variable = factor(case_when(variable == "Treating members of all religious groups equally" ~ "Treating members of\nall religious groups\nequally",
                                     variable == "More permissive policies toward undocumented immigrants" ~ "More permissive\npolicies toward\nundocumented immigrants",
                                     variable == "Preventing police use of force against peaceful protesters" ~ "Preventing police\nuse of force against\npeaceful protesters",
                                     variable == "Protection of media from government censorship" ~ "Protection of media\nfrom government\ncensorship",
                                     variable == "Affirmative action in university admissions" ~ "Affirmative action in\nuniversity admissions"),
                           levels = c("Treating members of\nall religious groups\nequally",
                                      "More permissive\npolicies toward\nundocumented immigrants",
                                      "Preventing police\nuse of force against\npeaceful protesters",
                                      "Protection of media\nfrom government\ncensorship",
                                      "Affirmative action in\nuniversity admissions"))) %>%
  # plot
  {ggplot(., aes(x = fct_rev(context), y = prop)) +
      geom_line(linewidth = 2, color = "gray") +
      geom_point(aes(color = religpew), size = 4) +
      geom_text(aes(label = percent(prop, accuracy = 1)), nudge_y = ifelse(.[["greater"]], 0.05, -0.05)) +
      facet_grid(rows = vars(variable), switch = "y") +
      ggtitle("Percent of respondents who support...") +
      coord_flip() +
      theme_bw() +
      theme(plot.title = element_text(size = 14, hjust = 0),
            plot.title.position = "plot",
            strip.placement = "outside",
            strip.background = element_blank(),
            strip.text.y.left = element_text(size = 12, face = "bold", angle = 0),
            legend.position = "top",
            legend.title = element_blank(),
            legend.text = element_text(size = 11),
            axis.title = element_blank(),
            axis.text = element_text(size = 11),
            axis.text.x = element_blank(),
            axis.ticks.x = element_blank(),
            panel.grid.minor = element_blank())}

# appendix

##### table B2 #####
main_ols <- lm(formula = diff2 ~ religion_ + education_ + immigrant_ + state_ + income_ + gender_ + age_ + biden2020_, data = iaas_wtd, weights = iaas_wtd[['variables']][['weight']])
summary(main_ols)

stargazer(main_ols, type = "text", title = "Linear Model (OLS)",
          covariate.labels = c("Religion: Hindu", "Religion: Other", "Education: Undergrad", "Education: Post-grad", "Immigrant", "Blue state", "Income over 100k", "Male", "Age 30 to 50", 
                               "Age over 50", "Voted for Biden in 2020"), 
          dep.var.labels = "Context Difference", star.cutoffs = c(0.05, 0.01, 0.001))

##### table B3 #####

iaas_wtd <- iaas_wtd %>%
  mutate(
    # demographic variables
    religion_ = factor(case_when(religpew == "Hindu" ~ "Hindu",
                                 religpew == "Muslim" ~ "Muslim",
                                 TRUE ~ "Other"), 
                       levels = c("Muslim", "Hindu", "Other")),
    education_ = factor(case_when(!str_detect(educ, "4-year|Post-grad") ~ "<Bachelor's",
                                  TRUE ~ educ)),
    state_lean = fec_states$lean[match(inputstate, fec_states$State)],
    state_ = factor(state_lean, levels = c("Red state", "Blue state")),
    income_ = factor(case_when(str_detect(faminc_new, "[0-9]{3},[0-9]{3} ") ~ ">$100,000",
                               faminc_new == "Prefer not to say" ~ NA_character_,
                               TRUE ~ "<$100,000")),
    gender_ = factor(case_when(gender == "Male" ~ "Male",
                               TRUE ~ "Not Male"), 
                     levels = c("Not Male", "Male")),
    age_ = factor(case_when(2020 - birthyr < 30 ~ "<30",
                            2020 - birthyr < 50 ~ "30-50",
                            TRUE ~ ">50"),
                  levels = c("<30", "30-50", ">50")), 
    biden2020_ = factor(case_when(B2_Q4 == "Joe Biden" ~ "Voted for Biden", 
                                  TRUE ~ "Did not vote for Biden"), 
                        levels = c("Did not vote for Biden", "Voted for Biden")))

# individual policy diff
iaas_wtd <- iaas_wtd %>%
  mutate(across(c(B4US_Q18:B4US_Q19, B4US_Q21:B4US_Q22, B5_Q11:B5_Q13, B5_Q15), list(code = ~case_when(.x == "Strongly support" ~ 1, 
                                                                                                       .x == "Somewhat support" ~ 2,
                                                                                                       .x == "Somewhat oppose" ~ 3,
                                                                                                       .x == "Strongly oppose" ~ 4))),
         across(c(B4US_Q20, B5_Q14), list(code = ~case_when(.x == "Strongly oppose" ~ 1,
                                                            .x == "Somewhat oppose" ~ 2,
                                                            .x == "Somewhat support" ~ 3,
                                                            .x == "Strongly support" ~ 4))),
         across(c(B5_Q11_code:B5_Q15_code, B5_Q14_code), ~-.x)) %>%
  mutate(press_diff = rowSums(select(.[['variables']], c(B4US_Q18_code, B5_Q11_code))), 
         mig_diff = rowSums(select(.[['variables']], c(B4US_Q19_code, B5_Q12_code))), 
         muslim_diff = rowSums(select(.[['variables']], c(B4US_Q21_code, B5_Q13_code))), 
         affirm_diff = rowSums(select(.[['variables']], c(B4US_Q20_code, B5_Q14_code))), 
         police_diff = rowSums(select(.[['variables']], c(B4US_Q22_code, B5_Q15_code))))

# run 5 models
redef_press <- svyglm(press_diff ~ religion_ + education_ + immigrant_ + state_ + income_ + gender_ + age_ + biden2020_, design=iaas_wtd)
redef_mig <- svyglm(mig_diff ~ religion_ + education_ + immigrant_ + state_ + income_ + gender_ + age_ + biden2020_, design=iaas_wtd)
redef_muslim <- svyglm(muslim_diff ~ religion_ + education_ + immigrant_ + state_ + income_ + gender_ + age_ + biden2020_, design=iaas_wtd)
redef_affirm <- svyglm(affirm_diff ~ religion_ + education_ + immigrant_ + state_ + income_ + gender_ + age_ + biden2020_, design=iaas_wtd)
redef_police <- svyglm(police_diff ~ religion_ + education_ + immigrant_ + state_ + income_ + gender_ + age_ + biden2020_, design=iaas_wtd)

stargazer(redef_press, redef_mig, redef_muslim, redef_affirm, redef_police, out = "Individual Policy Positions as DV.tex", title = "Individual Policies as DV", 
          covariate.labels = c("Religion: Hindu", "Religion: Other", "Education: Undergrad", "Education: Post-grad", "Immigrant", "Blue state", "Income over 100k", "Male", "Age 30 to 50", 
                               "Age over 50", "Biden voter (2020)"), 
          dep.var.labels = c("Media Censorship", "Immigration", "Religious Equality", "Affirm. Action", "Police Violence"), 
          star.cutoffs = c(0.05, 0.01, 0.001), type="text")

##### table C4 #####
iaas_wtd <- iaas_wtd %>%
  mutate(
    discr_ = factor(case_when(B4Filler_Q17_6 == "not selected" ~ "Have experienced discrimination", 
                              TRUE ~ "Have not experienced discrimination"), 
                    levels = c("Have not experienced discrimination", "Have experienced discrimination")))

# creating DV 
iaas_wtd <- iaas_wtd %>%
  mutate(across(c(B4US_Q18:B4US_Q19, B4US_Q21:B4US_Q22, B5_Q11:B5_Q13, B5_Q15), list(code = ~case_when(.x == "Strongly support" ~ 1, 
                                                                                                       .x == "Somewhat support" ~ 2,
                                                                                                       .x == "Somewhat oppose" ~ 3,
                                                                                                       .x == "Strongly oppose" ~ 4))),
         across(c(B4US_Q20, B5_Q14), list(code = ~case_when(.x == "Strongly oppose" ~ 1,
                                                            .x == "Somewhat oppose" ~ 2,
                                                            .x == "Somewhat support" ~ 3,
                                                            .x == "Strongly support" ~ 4))),
         across(c(B5_Q11_code:B5_Q15_code, B5_Q14_code), ~-.x)) %>%
  mutate(diff = rowSums(select(.[["variables"]], c(B4US_Q19_code, B4US_Q21_code, B5_Q12_code, B5_Q13_code))),
         diff2 = rowSums(select(.[["variables"]], ends_with("_code"))))

redef_discr <- svyglm(diff2 ~ religion_ + education_ + immigrant_ + state_ + income_ + gender_ + age_ + biden2020_ + discr_, 
                      design=iaas_wtd)
redef_discr %>% summary()

stargazer(redef_discr, type = "text", title = "Linear Model, Incl. Discrimination", 
          covariate.labels = c("Religion: Hindu", "Religion: Other", "Education: Undergrad", "Education: Post-grad", "Immigrant", "Blue state", "Income over 100k", "Male", "Age 30 to 50", 
                               "Age over 50", "Biden voter (2020)", "Have exp. discrimination"), 
          dep.var.labels = "Context Difference", star.cutoffs = c(0.05, 0.01, 0.001))

##### table D5 #####
ideo_mods <- svyglm(diff2 ~ religion_ + education_ + immigrant_ + state_ + income_ + gender_ + age_ + biden2020_, 
                    design = iaas_wtd, subset = iaas_wtd$variables$ideo7_client == "Moderate; middle of the road")
summary(ideo_mods)

stargazer(ideo_mods, type = "text", title = "Linear Model, Restricted to Moderates (ideo7)",
          covariate.labels = c("Religion: Hindu", "Religion: Other", "Education: Undergrad", "Education: Post-grad", "Immigrant", "Blue state", "Income over 100k", "Male", "Age 30 to 50", 
                               "Age over 50", "Voted for Biden in 2020"), 
          dep.var.labels = "Context Difference", star.cutoffs = c(0.05, 0.01, 0.001))

##### table E6 #####
iaas_wtd <- iaas_wtd %>% 
  mutate(
    # baseline/demographic variables
    religion_ = factor(case_when(religpew == "Hindu" ~ "Hindu",
                                 religpew == "Muslim" ~ "Muslim",
                                 TRUE ~ "Other"), 
                       levels = c("Muslim", "Hindu", "Other")),
    education_ = factor(case_when(!str_detect(educ, "4-year|Post-grad") ~ "<Bachelor's",
                                  TRUE ~ educ)),
    state_lean = fec_states$lean[match(inputstate, fec_states$State)],
    state_ = factor(state_lean, levels = c("Red state", "Blue state")),
    income_ = factor(case_when(str_detect(faminc_new, "[0-9]{3},[0-9]{3} ") ~ ">$100,000",
                               faminc_new == "Prefer not to say" ~ NA_character_,
                               TRUE ~ "<$100,000")),
    gender_ = factor(case_when(gender == "Male" ~ "Male",
                               TRUE ~ "Not Male"), 
                     levels = c("Not Male", "Male")),
    age_ = factor(case_when(2020 - birthyr < 30 ~ "<30",
                            2020 - birthyr < 50 ~ "30-50",
                            TRUE ~ ">50"),
                  levels = c("<30", "30-50", ">50")), 
    biden2020_ = factor(case_when(B2_Q4 == "Joe Biden" ~ "Voted for Biden", 
                                  TRUE ~ "Did not vote for Biden"), 
                        levels = c("Did not vote for Biden", "Voted for Biden")), 
    
    # engagement with india 
    infam_ = factor(case_when(B1_Q16 == "Yes" ~ "Family in India", 
                              TRUE ~ "No Family in India"), 
                    levels = c("No Family in India", "Family in India")), 
    engagement_ = rowSums(select(.[["variables"]], B4Filler_Q4_1:B4Filler_Q4_4) %>%
                            mutate(across(everything(), ~if_else(.x == "selected", 1, 0)))),
    # engagement_ := factor(case_when(engagement_ == 0 ~ "Low engagement",
    #                                 engagement_ == 1 ~ "Medium engagement",
    #                                 TRUE ~ "High engagement"),
    #                       levels = c("Low engagement", "Medium engagement", "High engagement")),
    # identity_ = factor(case_when(str_detect(B4Filler_Q8, "unimportant") ~ "Unimportant",
    #                              TRUE ~ B4Filler_Q8),
    #                    levels = c("Unimportant", "Somewhat important", "Very important")), 
    identity_ = case_when(B4Filler_Q8 == "Very unimportant" ~ 1, 
                          B4Filler_Q8 == "Somewhat unimportant" ~ 2, 
                          B4Filler_Q8 == "Somewhat important" ~ 3, 
                          B4Filler_Q8 == "Very important" ~ 4), 
    inpoleng_ = factor(case_when(B5_Q1 == "Most of the time" | B5_Q1 == "Some of the time" ~ "Politically engaged in India", 
                                 TRUE ~ "Not politically engaged in India"), 
                       levels = c("Not politically engaged in India", "Politically engaged in India")))

# B4Filler_Q4_1 Engaged in these activities: Traveled to India in the last one year
# B4Filler_Q4_2 Engaged in these activities: Communicate with friends and family in India at least once a month
# B4Filler_Q4_3 Engaged in these activities: Supported non-profit/developmental/philanthropic organizations in India in the last one year
# B4Filler_Q4_4 Engaged in these activities: Supported religious organizations in India in the last one year

# creating DV 
iaas_wtd <- iaas_wtd %>%
  mutate(across(c(B4US_Q18:B4US_Q19, B4US_Q21:B4US_Q22, B5_Q11:B5_Q13, B5_Q15), list(code = ~case_when(.x == "Strongly support" ~ 1, 
                                                                                                       .x == "Somewhat support" ~ 2,
                                                                                                       .x == "Somewhat oppose" ~ 3,
                                                                                                       .x == "Strongly oppose" ~ 4))),
         across(c(B4US_Q20, B5_Q14), list(code = ~case_when(.x == "Strongly oppose" ~ 1,
                                                            .x == "Somewhat oppose" ~ 2,
                                                            .x == "Somewhat support" ~ 3,
                                                            .x == "Strongly support" ~ 4))),
         across(c(B5_Q11_code:B5_Q15_code, B5_Q14_code), ~-.x)) %>%
  mutate(diff = rowSums(select(.[["variables"]], c(B4US_Q19_code, B4US_Q21_code, B5_Q12_code, B5_Q13_code))),
         diff2 = rowSums(select(.[["variables"]], ends_with("_code"))))

# a) baseline model with "family in india" var as control
lm_fam <- svyglm(diff2 ~ religion_ + education_ + immigrant_ + state_ + income_ + gender_ + age_ + biden2020_ + infam_, 
                 design=iaas_wtd)
lm_fam %>% summary()

stargazer(lm_fam, title = "Linear Model, Family in India", 
          covariate.labels = c("Religion: Hindu", "Religion: Other", "Education: Undergrad", "Education: Post-grad", "Immigrant", "Blue state", "Income over 100k", "Male", "Age 30 to 50", 
                               "Age over 50", "Biden voter (2020)", "Family in India"), 
          dep.var.labels = "Context Difference", star.cutoffs = c(0.05, 0.01, 0.001))


# b) baseline model with "engagement with india" var coded as CONTINUOUS
lm_engage <- svyglm(diff2 ~ religion_ + education_ + immigrant_ + state_ + income_ + gender_ + age_ + biden2020_ + engagement_, 
                    design=iaas_wtd)
lm_engage %>% summary()

stargazer(lm_engage, title = "Linear Model, Engagement with India", out = "LM (Engagement with India).tex", 
          covariate.labels = c("Religion: Hindu", "Religion: Other", "Education: Undergrad", "Education: Post-grad", "Immigrant", "Blue state", "Income over 100k", "Male", "Age 30 to 50", 
                               "Age over 50", "Biden voter (2020)", "Engagement with India"), 
          dep.var.labels = "Context Difference", star.cutoffs = c(0.05, 0.01, 0.001))


# c) baseline model with "importance of indian identity" var coded as CONTINUOUS
lm_id <- svyglm(diff2 ~ religion_ + education_ + immigrant_ + state_ + income_ + gender_ + age_ + biden2020_ + identity_, 
                design=iaas_wtd)
lm_id %>% summary()

stargazer(lm_id, title = "Linear Model, Importance of Indian Identity", out = "LM (Importance of Indian Identity).tex", 
          covariate.labels = c("Religion: Hindu", "Religion: Other", "Education: Undergrad", "Education: Post-grad", "Immigrant", "Blue state", "Income over 100k", "Male", "Age 30 to 50", 
                               "Age over 50", "Biden voter (2020)", "Indian Identity Important"), 
          dep.var.labels = "Context Difference", star.cutoffs = c(0.05, 0.01, 0.001))


# d) baseline model with "involvement in indian politics" var as is
lm_poleng <- svyglm(diff2 ~ religion_ + education_ + immigrant_ + state_ + income_ + gender_ + age_ + biden2020_ + inpoleng_, 
                    design=iaas_wtd)
lm_poleng %>% summary()

stargazer(lm_poleng, title = "Linear Model, Involvement in Indian Politics", 
          covariate.labels = c("Religion: Hindu", "Religion: Other", "Education: Undergrad", "Education: Post-grad", "Immigrant", "Blue state", "Income over 100k", "Male", "Age 30 to 50", 
                               "Age over 50", "Biden voter (2020)", "Politically Engaged in India"), 
          dep.var.labels = "Context Difference", star.cutoffs = c(0.05, 0.01, 0.001))

# combined model
stargazer(lm_fam, lm_engage, lm_id, lm_poleng, star.cutoffs = c(0.05, 0.01, 0.001))

##### table F7 #####
bluestates <- svyglm(diff2 ~ religion_ + education_ + immigrant_ + state_ + income_ + gender_ + age_ + biden2020_, 
                     design = iaas_wtd, subset = inputstate %in% c("Texas", "New Jersey", "New York", "California", "Florida", "Illinois"))
bluestates %>% summary()

stargazer(bluestates, type = "text", title = "Linear Model, Restricted to TX, NJ, NY, CA, FL, IL", 
          covariate.labels = c("Religion: Hindu", "Religion: Other", "Education: Undergrad", "Education: Post-grad", "Immigrant", "Blue state", "Income over 100k", "Male", "Age 30 to 50", 
                               "Age over 50", "Voted for Biden in 2020"), 
          dep.var.labels = "Context Difference", star.cutoffs = c(0.05, 0.01, 0.001))

##### table G8 #####

# with interaction subset to hindus and muslims only
order2 <- svyglm(diff2 ~ Mod_Block4And5_order + religion_ + education_ + immigrant_ + state_ + income_ + gender_ + age_ + biden2020_ +
                   Mod_Block4And5_order * religion_, design=iaas_wtd, subset = religion_ %in% c("Hindu", "Muslim"))
summary(order2)
stargazer(order2,  dep.var.labels = "Context Difference", star.cutoffs = c(0.05, 0.01, 0.001), type="text")

##### table H9 #####

iaas_wtd <- iaas_wtd %>% 
  mutate(
    arrivalyear = as.numeric(B1_Q5), 
    birthyear = as.numeric(birthyr), 
    yearsindia = arrivalyear - birthyear
  )

yearsin <- svyglm(diff2 ~ religion_*yearsindia + education_ + state_ + income_ + gender_ + age_ + biden2020_, 
                  design=iaas_wtd)
yearsin %>% summary()

stargazer(yearsin, type = "text", title = "Linear Model, Hindu x Years Spent in India", omit = "religion_Other:yearsindia", 
          covariate.labels = c("Religion: Hindu", "Religion: Other", "Years Spent in India (Year of Arrival - Birth Year)", "Education: Undergrad", "Education: Post-grad", "Blue state", "Income over 100k", 
                               "Male", "Age 30 to 50", "Age over 50", "Biden voter in 2020", "Religion: Hindu x Years Spent in India"), 
          dep.var.labels = "Context Difference", star.cutoffs = c(0.05, 0.01, 0.001))

##### table H10 #####

hindu_mig <- svyglm(diff2 ~ religion_*immigrant_ + education_ + state_ + income_ + gender_ + age_ + biden2020_, 
                    design=iaas_wtd)
hindu_mig %>% summary()


stargazer(hindu_mig, type = "text", title = "Linear Model, Hindu x Immigrant", omit = "religion_Other:immigrant_Immigrant", 
          covariate.labels = c("Religion: Hindu", "Religion: Other", "Immigrant", "Education: Undergrad", "Education: Post-grad", "Blue state", "Income over 100k", 
                               "Male", "Age 30 to 50", "Age over 50", "Biden voter in 2020", "Religion: Hindu x Immigrant"), 
          dep.var.labels = "Context Difference", star.cutoffs = c(0.05, 0.01, 0.001))

##### table I11 and I12 ######

iaas$order_USfirst <- NA
iaas$order_USfirst[iaas$Mod_Block4And5_order=="Normal order: Block4US, Block4Filler, then Block5"] <- 1
iaas$order_USfirst[iaas$Mod_Block4And5_order=="Reverse order: Block5, Block4Filler, then Block4US"] <- 0
table(iaas$order_USfirst)
iaas_wtd <- update(iaas_wtd, order_USfirst = iaas$order_USfirst)


levels_liberal <- c("Strongly support", "Somewhat support", "Somewhat oppose", "Strongly oppose")

# Recode function
recode_liberal <- function(x, reverse = FALSE) {
  x <- factor(x, levels = levels_liberal)
  if (reverse) {
    return(as.numeric(x))          # conservative-to-liberal coding
  } else {
    return(5 - as.numeric(x))      # liberal-to-conservative flipped
  }
}

# Apply to variables
iaas_wtd <- iaas_wtd %>%
  mutate(
    # U.S. items
    Q18_lib = recode_liberal(B4US_Q18, reverse = TRUE),   # press credentials
    Q19_lib = recode_liberal(B4US_Q19, reverse = TRUE),   # ICE raids
    Q20_lib = recode_liberal(B4US_Q20, reverse = FALSE),  # affirmative action (race)
    Q21_lib = recode_liberal(B4US_Q21, reverse = TRUE),   # Muslim travel ban
    Q22_lib = recode_liberal(B4US_Q22, reverse = TRUE),   # force on BLM protesters
    
    # India items
    Q11_lib = recode_liberal(B5_Q11, reverse = TRUE),     # sedition laws vs reporters
    Q12_lib = recode_liberal(B5_Q12, reverse = TRUE),     # NRC
    Q13_lib = recode_liberal(B5_Q13, reverse = TRUE),     # CAA
    Q14_lib = recode_liberal(B5_Q14, reverse = FALSE),    # caste-based admissions
    Q15_lib = recode_liberal(B5_Q15, reverse = TRUE)      # force on anti-CAA protesters
  )

iaas <- iaas %>%
  mutate(
    # U.S. items
    Q18_lib = recode_liberal(B4US_Q18, reverse = TRUE),   # press credentials
    Q19_lib = recode_liberal(B4US_Q19, reverse = TRUE),   # ICE raids
    Q20_lib = recode_liberal(B4US_Q20, reverse = FALSE),  # affirmative action (race)
    Q21_lib = recode_liberal(B4US_Q21, reverse = TRUE),   # Muslim travel ban
    Q22_lib = recode_liberal(B4US_Q22, reverse = TRUE),   # force on BLM protesters
    
    # India items
    Q11_lib = recode_liberal(B5_Q11, reverse = TRUE),     # sedition laws vs reporters
    Q12_lib = recode_liberal(B5_Q12, reverse = TRUE),     # NRC
    Q13_lib = recode_liberal(B5_Q13, reverse = TRUE),     # CAA
    Q14_lib = recode_liberal(B5_Q14, reverse = FALSE),    # caste-based admissions
    Q15_lib = recode_liberal(B5_Q15, reverse = TRUE)      # force on anti-CAA protesters
  )

# Create additive indices for U.S. and India 
iaas_wtd <- iaas_wtd %>%
  mutate(
    us_liberal = Q18_lib + Q19_lib + Q20_lib + Q21_lib + Q22_lib,
    india_liberal = Q11_lib + Q12_lib + Q13_lib + Q14_lib + Q15_lib
  )

iaas <- iaas %>%
  mutate(
    us_liberal = Q18_lib + Q19_lib + Q20_lib + Q21_lib + Q22_lib,
    india_liberal = Q11_lib + Q12_lib + Q13_lib + Q14_lib + Q15_lib
  )

#liberalism score for block they saw first 
iaas_wtd <- iaas_wtd %>%
  mutate(
    liberalism_first = case_when(
      order_USfirst == 1 ~ us_liberal,   # saw US first → use US score
      order_USfirst == 0 ~ india_liberal,# saw India first → use India score
      TRUE ~ NA_real_
    )
  )

# table I 11

order_new <- svyglm(liberalism_first ~ order_USfirst + religion_ + education_ + immigrant_ + state_ + income_ + gender_ + age_ + biden2020_, 
                    design=iaas_wtd)
summary(order_new)
stargazer(order_new, star.cutoffs = c(0.05, 0.01, 0.001), type="text")

# table I 12
order_new2 <- svyglm(liberalism_first ~ order_USfirst * religion_ + education_ + immigrant_ + state_ + income_ + gender_ + age_ + biden2020_, 
                     design=iaas_wtd)
summary(order_new2)
stargazer(order_new2, type="text", star.cutoffs = c(0.05, 0.01, 0.001))
